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The Fluctuation Relation (FR) is an asymptotic result on the distribution of certain observables 
averaged over time intervals r as t — » oo and it is a generalization of the fluctuation-dissipation 
theorem to far from equilibrium systems in a steady state which reduces to the usual Green-Kubo 
(GK) relation in the limit of small external non conservative forces. FR is a theorem for smooth 
uniformly hyperbolic systems, and it is assumed to be true in all dissipative "chaotic enough" 
systems in a steady state. In this paper we develop a theory of finite time corrections to FR, needed 
to compare the asymptotic prediction of FR with numerical observations, which necessarily involve 
fluctuations of observables averaged over finite time intervals r. We perform a numerical test of FR 
in two cases in which non Gaussian fiuctuations are observable while GK does not apply and we 
get a non trivial verification of FR that is independent of and different from linear response theory. 
Our results are compatible with the theory of finite time corrections to FR, while FR would be 
observably violated, well within the precision of our experiments, if such corrections were neglected. 
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I. INTRODUCTION 

Anosov systems and the fluctuation theorem - The fluctu- 
ation theorem concerns fluctuations of phase space con- 
traction in reversible hyperbolic (Anosov) systems. If 
time evolution is described by a differential equation 
on phase space M: x = X{x), x G M, or by a map 
S : X ^ S{x) of M one deflnes the phase space con- 
traction as, respectively, <j{x) — — divX(a;) or o{x) — 
— log I det i9a;<S'(a;)|. Reversibility means that there is a 
metric-preserving map / of M such that IS = S^^I if S 
is the time evolution over a certain time t {e.g. t — 1). 
// the system is Anosov, that is if M is compact and S 
is smooth and uniformly hyperbolic, see [1-5], the points 
X will have a well deflned SRB distribution fisrb, [5], i.e. 
almost all points w.r.t. the volume measure will evolve 
so that all smooth observables will have a well deflned 
average equal to the integral over the SRB distribution. 
Hence, in particular, the time average of the function 
cr(x) will be asymptotically given by the spatial average 
w.r.t. the SRB distribution. In the case of discrete time 
maps: 

1 f 

a+ = lim - a{S^{x)) = cr dfigrb ( cr )srb (1) 
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If a+ > 0, let: 

Analogous definitions are given in the continuous time 
case. The function p{x) will have average ( p )srb = 1 
and distribution TTridp) such that 

TTriip e A}) = e^™'''^''^^'^~(f)+°(^) , (3) 

where the correction at the exponent is o(r) w.r.t. t as 
T — > oo. The following Fluctuation Relation, discovered 
in a numerical simulation in [6] and formulated as a the- 
orem for Anosov systems in [2], holds: 

Coo [p] = Coo i-p) + pcr+ for all IpI < p* (4) 

where cx) > p* > 1 is a suitable (model dependent) con- 
stant that, in general, is different from the maximum over 
T and X oip{x); note also that Eq. 4 is (strictly speaking) 
meaningless in the equilibrium cases in which the system 
is Hamiltonian and reversibility is the usual velocity sign 
change: because of the division by cr+ = in Eq. 2. 

The chaotic hypothesis - Hyperbolicity is a paradigm 
for disordered systems similar to the small oscillations 
paradigm used for ordered motions: it does not hold ex- 
actly in essentially all the physically interesting systems. 
The Chaotic Hypothesis [1-3, 7, 8] is that nevertheless 
one can assume that chaotic motions (in the sense of 
motions with at least one positive Lyapunov exponent) 
exhibit some average properties of truly hyperbolic mo- 
tions. This hypothesis is a natural generalization of the 
ergodic hypothesis, i.e. of the assumption that systems 
of many particles at equilibrium are well described on 
average by the microcanonical (or by the Gibbs) distri- 
bution, even if they are not really (or they are not proven 
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to be) ergodic. A consequence of the Chaotic Hypothesis 
is that (dissipative) deterministic chaotic reversible mo- 
tions should have fluctuations of phase space contraction 
satisiying Eq. 4. 

One interesting example of such motions is given by 
a system of A'^ interacting particles in d dimensions sub- 
jected to nonconservative forces and kept in a station- 
ary state by a reversible mechanical thermostat. It will 
be defined by a differential equation x = Xe{x) where 
X = (9) q) S ii^''^ = M {phase space) and 



mq 



(5) 



where m is the mass of the particles, f{q) describes the 
internal (conservative) forces between the particles and 
9^{q) represents the nonconservative "external" force 
acting on the system. Finally, 9e{q,q) is a mechanical 
force that prevents the system from acquiring energy in- 
definitely: this is why we shall call it a mechanical ther- 
mostat. Systems belonging to this class are frequently 
used as microscopic models to describe nonequilibrium 
stationary states induced by the application of a driving 
force (temperature or velocity gradients, electric fields, 
etc.) on a fluid system in contact with a thermal bath, 
[8, 9]. In this context, the phase space contraction rate 
a{x) has been identified (setting fcs = 1) with the en- 
tropy production rate [1, 6-8], the variable p(a;) is defined 



as 



p{x) 



TO+ Jo 



dt a{Stx) 



(6) 



(where x{t) = Stx is the solution of Eq. 5 with initial 
datum x(0) — x) and the fluctuation relation has been 
successfully tested in several numerical simulations [6, 
10-15]. Having defined the notion of entropy production 
rate one can define a "duality" between fluxes J and 
forces using a{x) as a "Lagrangian" [7]: 



J{E,x) 



da{x) 
dE 



(7) 



In the limit E_^Q, i.e. close to equilibrium, the fluctua- 
tion relation leads to Onsager's reciprocity and to Green- 
Kubo's formulas for transport coeSicients [16, 17]: 



lim 



{Ji)l 

E, 



Jo 



dt {Mt)Jj{0))E=o (8) 



Gaussian distributions - If the distribution iTT-ip) is Gaus- 



sian. 



T^r{p) 



oc exp 



2(5? 



from the fluctuation re- 



lation one can derive an extension of the Green-Kubo 
relation, i.e. of Eq. 8, to finite forces. 

Indeed, the fluctuation relation for a Gaussian distri- 
bution implies that the dispersion 5^ of p around its 
average (equal to 1) is S'^ =2/(7+ which is, in such case, 
an extension of a Green-Kubo formula to non zero fields. 
One sees this by considering, for instance, cases in which 



(j{x) is linear in E (as it will be in the cases that we 
study numerically below). Using time- translation invari- 
ance one can show that 

2 f°° 

^lo = -T dt {{a{t) - a+)(a(0) - a+))E (9) 

K Jo 

and from the fluctuation relation (5^(7+ =2 

/•OO 

<7+ = / dt {{a{t) - a+){a{0) - a+))E (10) 
Jo 

Substituting a{t) — EJ{t) in the latter expression, one 
obtains the relation 



E 



dt [(J(i)J(0))B - {J)\ 



(11) 



valid, subject to the Gaussian assumption, also for E ^Q. 

The leading order in E of the latter relation [linear 
response) is the Green-Kubo formula for the equilibrium 
transport coefficient, Eq. 8. 

Numerical verification of the chaotic hypothesis - The 
simplest check of the applicability of the Chaotic Hy- 
pothesis is a check of the fiuctuation relation: of course 
even if the check has a positive result this will not "prove" 
the hypothesis but it will at least add confidence to it. A 
rather stringent test of the fluctuation relation would be 
a check which cannot be reduced to a kind of Green-Kubo 
relation; this requires at least one of the two following 
conditions to be satisfied: 

1. the distribution 'Krip) is distinguishable from a 

Gaussian, or 

2. deviations from the leading order in E in Eq. 11, 
i.e., deviations from the Green-Kubo relation, are 
observed. 

This is very hard to obtain in numerical simulations of 
Eq. 5 for the following reasons: 

1. to observe deviations from linearity in Eq. 11 one 
has to apply very large forces E, then (7+ is very 
large and it becomes very difficult to observe the 
negative values of p{x) that are needed to compute 
Coo(-p) in Eq. 4; 

2. deviations from Gaussianity in TTrip) are observed 
only for values of p that differ significantly (of the 
order of 25 oo) from 1 and, again, it is very difficult 
to observe such values of p. 

Due to the limited computational resources available in 
the past decade, all numerical computations that can be 
found in the literature on systems described by Eq. 5 
found that the measured distribution tt,- [p) could not be 
distinguished from a Gaussian distribution in the interval 
of p accessible to the numerical experiment [6, 10, 11, 15]. 

The purpose of the present paper is to test the fluc- 
tuation relation, in a numerical simulation of a system 
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described by Eq. 5, for large applied force when devia- 
tions from linearity can be observed, and the distribution 
TTr (p) is appreciably non-Gaussian. This has become pos- 
sible thanks to the fast increase of computational power 
in the last decade. However, it is still very difficult to 
reach values of r which can be confidently regarded as 
"close" to the asymptotic limit r — > c»; thus to interpret 
our results we develop a theory of the o(l) corrections to 
the function Coo(p) in order to extract the limiting func- 
tion (p) from the numerical data. Taking into account 
the latter finite time corrections, we successfully test the 
fluctuation relation for non-Gaussian distributions and 
beyond the linear response theory. 

The paper is organized as follows: section II is devoted 
to the theory of finite r corrections to the large deviation 
function (p) that is needed in the analysis of the data; 
in section III we present the model and the details of 
the numerical simulation; in section IV we present the 
details of the data analysis; finally, in section V and VI 
we report the result of our simulations. 



II. FINITE TIME CORRECTIONS TO THE 
FLUCTUATION RELATION 

In the present section we describe a strategy to study 
(in principle constructively) the 0(1) corrections in the 
exponent of Eq. 3. The theory we propose will hold as- 
suming that the time evolution is hyperbolic so that it can 
be applied to physical systems only if the chaotic hypoth- 
esis is accepted. For simplicity we consider only the case 
of discrete time evolution via a map S. 



SRB measure, symbolic dynamics and 
statistical mechanics 



We study the distribution of p at fixed r via its Laplace 
transform {characteristic function) Zr{X)- 



Zr{\) = — log(e 

T 



(12) 



The main consequence of the hyperbolicity is, [2-4, 17- 
19], that one can find a symbolic representation of the 

points of M in terms of sequences e = (eOi^-oo finitely 
many digits e = 1, . . . , fc subject only to a simple hard 
core restriction, namely Te^^e^_^j = 1 if T is a matrix 
(compatibility matrix) with entries or 1 and such that 
rj^, > for some iV > and all e, e' {mixing condition). 
Moreover in such a representation the dynamics becomes 
simply the left shift, i.e. if s{x) represents x then S{x) is 
represented by the sequence e shifted to the left by one 
unit. 

The key remarks are 

1. Smooth observables on phase space can be repre- 
sented by short range potentials: in the case of the 
observable a{x) this means that there are functions 



sx{£x) defined for all intervals X = {a, . . . , a+2n+ 
1) and Sx = {sai ■ ■ ■ ■,£a+2n+i)i translationally in- 
variant sx = sx+b and exponentially decaying on 
time scale {i.e. \sx{sx)\ < Ce"*^" for some 
C,K > 0), such that 

a{x) = s{e{x)), s{s) = (13) 

XoO 

where the sum is over the intervals X centered at 
the origin (noted by XoO). Another important 
smooth observable is the expansion rate L{x) de- 
fined as the logarithm of the determinant of the 
linearization matrix dS{x) {i.e. the Jacobian ma- 
trix of the map) restricted to the unstable manifold: 
L{x) = logdet9S'(a;)„. This is also expressible via 
an exponentially decaying potential $: 

L{x)=l{e{x)), m^Y.'^^i^x) (14) 

2. The SRB distribution, represented as a distribu- 
tion over the (compatible) symbolic sequences e, is 
a Gibbs state for the short range potential $ = 
{<^x{§.x)) defined in Eq. 14, i.e. 



srh 



lim 

ii->oo 



■E 



(15) 



where Kr = (— i?, . . . ,R) C Z, the sums extend 
over compatible configurations e = {s-r, . . . , er) 
{i.e. with Te^.,,^.^, = 1 for j = -R, . . . , R - 1), 
and F{e') is an arbitrary smooth observable de- 
fined on phase space regarded as a function on the 
symbolic sequences and evaluated at a sequence e' 
which extends (rather arbitrarily) s to an infinite 
compatible sequence by continuing e to the right 
with a sequence and to the left with a sequence 
into e' {e^ ,£,£>) so that e< depends only on 
the symbol E-r and £^ depends only on the symbol 
er: see [18, 20, 21]. 

The surprising reduction of the problem of studying the 
SRB distribution to that of a Gibbs distribution for a one 
dimensional chain with short range interaction (this is the 
physical interpretation of Eq. 15) generated the possibil- 
ity of studying more quantitatively at least some of the 
problems of nonequilibrium statistical mechanics outside 
the domain of "nonequilibrium thermodynamics", [22]. 



B. Finite time corrections to the characteristic 
function 

The characteristic function Zt{X) of p, see Eq. 12, can 
therefore be computed as 

g-ExcAH Exo [0,T-1] «x(£x) 

e-r..(A) ^ lij^ ^ ^ ^ 



(16) 
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This means that it is the (hmit as i? — > oo of the) ratio 
between the partition functions Zji{^) of a Gibbs dis- 
tribution in Afl with potential $ (the denominator) and 
the partition function Zu{^, As) with the same potential 
modified in the finite region [0, r — 1] C Z by the addition 
of a potential Xsxi^x)- 

The one dimensional systems are very well understood 
and the above is a well studied problem in statistical 
mechanics, known as a finite size corrections calculation. 
For instance it can be attacked by cluster expansion, [21]; 
this is a technique to deal with the average of the expo- 
nential of a spin Hamiltonian which is defined in terms of 
potentials (l)x exponentially decaying with rate k, such as 
those appearing in the numerator and in the denominator 
of Eq. 16. It allows us to represent them as: 



-Ex 



CAR 



(17) 



where (px are new effective potentials, explicitly com- 
putable in terms of suitable averages of products of 

cbxisxY^- and which can be proven to be still expo- 
nentially decaying with the diameter of X with a rate 

0<k' <K. 

In particular, a representation like Eq. 17 allows us 
to rewrite the partition function in the denominator of 
Eq. 16 as: 



Zr{^) = exp {2R + l)/oo($) - Coo($) + 0(e-'='«) 



and the one in the numerator as 



(18) 



Zr{^, As) = exp {2R+1- t)/oo($) + t/oo($ + As) 



- Coo($) - 5cx>(A) + 0(e-«'« + e-«'^) 



Therefore 



(19) 



^x(A) = /oo(*) - ,/oo(1> + As) + + 0{e-'^) 

T 



'M,^^X) + ^-^ + 0{e--'^) 



(20) 



The function Zaa{X) is convex in A and the functions 

.goc(A) and Zt{\) are analytic in A (a consequence of the 
1-dimensionality and of the short range nature of the 
SRB distribution): namely 500 (A) = X+^g'i^ )? + . . . 

and Zt{\) = z^^'X + ^z^^X"^ + . ■ ■ and the coefficients of 
their expansion in a power series of A can be expressed 
in terms of correlation functions of 17(0;). For instance, 
from Eq. 12 and using the translational invariance of the 



SRB measure, 

T-l 



.(1) 



.(2) 



3=0 



T-l 



srb 



3=0 

T-l 

■ E 

k=-T+l 



T 



3=0 



k=0 



(21) 



where {a{S''x)a{x))c = {a{S^x)a{x)) srb — o'+- Using 
Eq. 20, ,9cc(A) = limT-^oo T[zr(A) — -2:oo(A)], and the an- 
alyticity of 2^(A), we have = lim^^oo '''[•^t''^ — 
Since the connected correlation hmction {g{S'^x)cf{x))c 
decays exponentially for fc ^ 00, we obtain 



00 

9t^= E m'y{S^x)a{x)), 



(22) 



fc=— 00 



C. Finite time corrections to Coo(p) 

A direct measurement of Zr (A) from the numerical data 
is difficult. What is really accessible to numerical ob- 
servation are the quantities ^\ogiTT-{{p G A}) in Eq. 3 
because the measured values oip are used to build an his- 
togram obtained by dividing the p-axis into sufficiently 
small bins A and counting how many values of p fall in 
the various bins. Let us choose the size of the bins A as 
|A| = O^Et/t), with Er a small parameter which will be 
eventually chosen Er = o(l), see Appendix A for a dis- 
cussion of this point. Let also pa be the center of the bin 
A. An application of a local form of central limit theo- 
rem, discussed in Appendix A. shows that the following 
asymptotic representation of 7rx({p & A}) holds: 



7r^({pe A}) = e^^-(f^)(l + o(l)) 



(23) 



where Ct(pa) can be interpolated by an analytic function 
of p, satisfying the equation 



Ct(p) = -Zr{Xp) + Xppa+ 



2t 



log 



27r 



(- 



4'(Ap) 



(24) 

and Xp is the inverse of p(A) = z^(A)/cr+. 

Using the previous equations, we now compute the low- 
est order finite time correction to Coo (p) around the max- 
imum. 

We rewrite Cr(p) as Cr(p) = Coo(p) + ^ + 0{^). 
By the analyticity of Cr(p)) we can write Coo (p) , Too (p) 



around p = 1 in the form: Coo(p) = ^C&^ip ~ 1)^ 
hC^^Hp - 1)^ + . . . and 7oo(p) = 7^^ + 7^Vp - 1) + • 



+ 
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Up to terms of order (p — 1)^ and higher in the series 
for -fooip) we can rewrite: 

<t(p) = 

T T V / 

(25) 

Thus, the finite time corrections to Coo(p) around its max- 
imum begin with a shift of the maximum at 



Po 



/oo 



(2) 



(26) 



To apply the latter result we need to compute 7,^^ in 
terms of observable quantities. And, in order to compute 
7,^^ we apply Eq. 24. First of all, we note that Xp is 
determined by the condition 

pa+ = 4(Ap) = a+ + 4'(0)Ap + O(A^) (27) 

where we used Eq. 21 and Eq. 22. Then, Xp = + 
0{{p — 1)^). Substituting this result into Eq. 24 and 
equating the terms of order 0(^^) at both sides we find: 



^(1) 

loo 



1 (3) 



The last equation can also be rewritten as: 



v(l) 



43) 



2c; 



(2) 



(28) 



(29) 



This can be proven recalling that (^^^ and Coo'' are deriva- 
tives of Coo(p) in p = 1, that can be obtained by difi'eren- 
tiating w.r.t. A (two or three times, respectively) the def- 
inition Coo (2:^ (A)/(j+) = — Zoo(A) + Az^(A) and comput- 
ing the derivatives in A = recalling that 2;^(0)/cr+ = 1. 
Plugging Eq. 29 into Eq. 26 we finally get 



Po = 1 - 



^(3) 
S>oc 



2r(C^^) 



(30) 



that is the main result of this section. The key point is 
that the moments Coo'' and Coo'' in Eq. 30 are quantities 
that can be measured from our empirical data (within 
an 0{t~^) error). We then have a verifiable prediction 
on the expected shift of the maximum at finite r. Our 
data agree very well with this prediction, see Fig. 2 and 
corresponding discussion in sec. IV below. 
Substituting Eq. 30 in Eq. 25, we finally find: 



CM = Vrip) + + o(r-i) 



1)2. 



(31) 



where rjr{p) is defined as 

(0) / 

, .def 7oo , ^ / 
VAp) = = +Cr \P 



>(3) 



2r(Ci 



(32) 



The key point of the above discussion was the validity of 
Eq. 23-24; see Appendix A for their derivation. 



D. Remarks 

(1) The shift away from 1 of the maximum of the func- 
tion Ct{p) at finite r, expressed by the second term in 
Eq. 32, is due to the asymmetry of the distribution tTt{p) 
around the average value p — \] consequently, it is pro- 
portional, at leading order in t~^, to Coo^ which is indeed 
a measure of the asymmetry of Coo {p) around p = 1 . This 
shift would be absent in the case of a symmetric distri- 
bution {e.g., a Gaussian) and for this reason it was not 
observed in previous experiments [6, 10, 11, 15]. 

(2) The error term in the r.h.s. of Eq. 23 is o(l) w.r.t. 
T and it does not affect the computation of 700 (p)- It is 
then clear that with a calculation similar to that we per- 
formed, one can get equations for the coefficients 0{X'^) 
in the exponents of Eq. 23; in this way one can iteratively 
construct the whole sequence of coefficients 7,^' defining 
the power series expansion of 7oo(p)- 

(3) In models with continuous time evolution the quan- 
tity (7-1- is not dimensionless but it has dimensions of in- 
verse time: in such cases one can imagine that one is 
still studying a map which maps a system configuration 
at a time when some prefixed event happens in the sys- 
tem (typically a "collision") into the next one in which 
a similar event takes place. If tq is the average time in- 
terval between such events then rocr+ will play the role 
played by 0"+ in the; (lisc;rete time case: it will be the adi- 
mensional parameter entering the estimates of the error 
terms. 

Note that the coefficients g^a' are of order and 
their size is necessarily estimated by (the adimensional) 
entropy production to the fc-th power. Then, in the con- 
tinuous time case, the choice of tq affects the estimates 
of the remainders, because it affects the size of the adi- 
mensional parameter Tocr+; and the size of the mixing 
time (that is connected with the estimated range of de- 
cay of the potentials, see [21]). The natural (and physi- 
cal) choice for tq is the mixing time. Consistently with 
this remark, at the moment of constructing numerically 
the distribution function for the entropy production rate 
averaged over a time r, we will always consider time in- 
tervals of the form r = tqu, n > 1, see section IIIC 
below. 



III. THE MODEL 

We consider a system of TV classical particles of equal 
mass m in dimension rf; they are described by their po- 
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sition qi and momenta pi = mqi, {pi,qi) G R^'^, i = 
1, . . . , N. The particles are confined in a cubic box of 
side L with periodic boundary conditions. Each particle 
is subject to a conservative force, fi{q) = —dq-V{q), and 
to a nonconservative force Ei that docs not depend on 
the phase space variables. The force is locally con- 
servative but not globally such due to periodic bound- 
ary conditions. The mechanical thermostat is a Gaussian 
thermostat [9], Oi{p,q) = —Q:{p,q)pi, and the function 
a{p, q) is defined by the condition that the total kinetic 

energy Kip) = ^kbl^ = 2k should be a constant 

{isokinetic ensemble). The equations of motion are: 



q^ = ^ ^ 

Pi = Mq) + Ei - a{p,q) Pi , 



Prom the constraint ^ = one obtains 



a{p, q) 



(33) 



(34) 



A. Entropy production rate 

The total phase space volume contraction rate for this 
system is given by: 



= dNa{p,q) + -^pi = {dN - 1) a{p,q) . 



(35) 



Defining the kinetic temperature, T = 2K{p)/{dN — 1), 
[9] , the phase space contraction rate can be rewritten as 



(^{P, q) 



EiEiqi-V 



(36) 



The first term is the power dissipated by the external 

force divided by the kinetic temperature, and can be 
identified with the entropy production rate [6, 7, 9]. The 
second term is the total derivative w.r.t. time of the 
potential energy divided by the temperature: this term 
does not affect the validity of the Fluctuation Relation 
in the asymptotic limit t ^ 00, as total derivatives give 
a contribution 0{t~^) in p{x) [7, 23], hence they do not 
contribute to (00 ; however it has effect on the distribution 
of fluctuations over a finite time r and its influence on 
the numerical computations has been recently discussed 
in detail [15]. The most convenient thing to do, in order 
to have a finite time distribution that approximates in 
the best possible way the asymptotic distribution of fluc- 
tuations, is to study the distribution of fluctuations for 
the entropy production rate s, where s is identified with 
a minus the total derivative term —V/T in Eq. 36: 



s{p, q) = 



Ei Ei qi 



(37) 



From now on we will call C,oo{p) and Ct(p) the distribu- 
tions for the fluctuations of the entropy production rate 
s averaged over infinite or finite time, respectively. These 
will be the objects we will measure and use from now on. 

In order to define the current J{x,E), let us rewrite 
Ei = Eui, where Ui is a (constant) unit vector that spec- 
ifies the direction of the force acting on the i-th particle. 
Then, according to Eq. 7, 



J {p, q) 



da_ 
dE 



T,i qi 



(38) 



B. Discretization of the equations of motion 

To perform the numerical simulation, one has to write 
the equations of motion in a discrete form. One possi- 
bility is to use the Verlet algorithm [24]; for Hamiltonian 
equations of motion {i.e., E = and a = 0) 



J ~ m ' 

\ Pi = fi{q) , 

the Verlet discretization has the form 



qiit + dt) 
Pi{t + dt) 



q,{t) + sMdt^ 

P^it) + ^[f^{t) 



\f^{t)dt^ , 
-f,{t + dt)]dt , 



(39) 



(40) 



where dt is the time step size. This discretization ensures 
that the error is 0{dt'^) on the positions qi{t) in a single 
time step. The implementation of this algorithm on a 
computer is discussed in detail in [24]. 

However, this method requires the forces fi{t) to de- 
pend only on the positions and not on the velocities: 
hence, it has to be adapted to Eq.s 33. This has been 
done in the following way. We write the discretized equa- 
tions as 




I^dt 

rn 

E, - a{t)p,{t)\dt^ 



= q^{t) 

+h[Mt) 

= Pi{t) + Ei + ^[fi{t) + fi{t + dt) 
-a{t)pi{t) -a{t + dt)pi{t + dt)] dt , 



(41) 



with the same error as in the standard Verlet discretiza- 
tion. We store in the computer, at time t, the positions 
qi{t^, the momenta Pi (f), the forces fi{t?), and the Gaus- 
sian multiplier a{t). Then, we perform the following op- 
erations: 

1. we calculate the new positions qi{t + dt) using the 
first equation; 

2. using the new positions we calculate the new forces 
fi {t + dt) (the conservative forces depend only on 
the positions); 

3. we calculate the quantity = pi{t) + Ei + i [fi{t) + 
fi{t+dt)—a{t)pi{t)] dt and we observe that pi{t+dt) 
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can be expressed in terms of the (known) and the 
(unknown) a{t + dt) as 



Pi{t + dt) 



a{t + dt)dt/2 



(42) 



4. substituting Eq. 42 in the definition of a{t + dt), 
Eq. 34, we get a self-consistency equation for a{t + 
dt), whose solution is 



a{t + dt) 



ceo 



ao 



1 - aodt/2 



■ dt) 



(43) 



5. substituting Eq. 43 in Eq. 42 we calculate pi(t + (it). 

This procedure allows us to calculate the new positions, 
momenta, forces, and a, at time t + dt according to 
Eq.s 41 without approximations, defining a map S such 
that {p{t + dt),q{t + dt)) = S{p{t),qit)). 

Our (discrete) dynamical system will be defined by the 
map S{p, q) and will approximate the differential equa- 
tions of motion, Eq. 33, with error 0{dt'^) for the posi- 
tions and 0{dt^) for the velocities. 
The map S satisfies the following properties: 



1. it is reversible, i.e. 
defined hy I{p,q) = 



it exists a map I{p, q) (simply 
: {-p,q)) such that"75 = S''^!; 



2. in the Hamiltonian case {E = and a = 0, Eq.s 39) 
it is volume preserving. 

The first property ensures that assuming the Chaotic Hy- 
pothesis the Fluctuation Relation holds for the map S. 
The second property ensures that at equilibrium the dis- 
cretization algorithm conserves the phase space volume. 



C. Details of the simulation 

In the simulation, we chose the external force of the 
form Ei = Eui, where the unit vectors Ui were parallel 
to the X direction but with different orientation: half of 
them were oriented in the positive direction, and half in 
the negative direction, i.e. Ui = (— l)*x, in order to keep 
the center of mass fixed. We considered two different 
systems, selecting interaction potentials widely used in 
numerical simulations (for the purpose of making easier 
possible future independent checks and rederivations of 
our results): 

1. {model I) the first investigated system is made by 
N = 8 particles of equal mass m in d = 2. The 
interaction potential is a sum of pair interactions, 
— Si<j ~ P^^^ interaction 

is represented by a WCA potential, i.e. a Lennard- 
Jones potential truncated at the minimum: 



v{r) 




7)" -iff 



+ e, r<^a ; 
r > . 



The reduced density was p = Na"^ /L"^ = 0.95 (that 
determines L), the kinetic temperature was fixed 
to T = 4e and the time step to dt = O.OOlto, where 
to = yma^Je. In the following, all the quantities 
will be reported in units of m, e and cr (LJ units). 
This system was already studied in the literature, 
see e.g. [6, 25]. We investigated different values of 
the external force E ranging from E = Qto E = 25. 

2. [model II) the second system is a binary mixtme of 
N=20 particles (16 of type A and 4 of type B), of 
equal mass m, in d = 3, interacting via the same 
WCA potential of model I; the pair potential is 



"a/3 



(r) 



0, 



6 

+ ea/3 , 

r > \plact0 



a and /? are indexes that specify the particle species 
(a,/? G {A,B}). The parameters entering the po- 
tential are the following: oab = O.Sctaa; (^bb = 
O.SSctaa; cab = I-Scaa; ess = 0.5e^^. Similar 
potentials have been studied, [26, 27], as models 
for liquids in the supercooled regime (i.e., below 
the melting temperature). For this system the LJ 
units are to, caa, and aAA', the miit of time is 
then to = \/ TTicr'^j^/eAA- The reduced density was 
p = Ncr\^/L^ = 1.2 and the integration step was 
dt = O.OOlto- The unit vectors are chosen such 
that half of the A particles and half of the B par- 
ticles have positive force in the x direction, and 
the remaining particles have negative force in the x 
direction. For this system we investigated different 
values of external force E e [0, 10] and temperature 
Te [0.5,3]. 

For each system and for each chosen value of T and E, 
we simulated a very long trajectory (~ 2 • lO^di) start- 
ing from a random initial data; we recall that in both 
systems we chose dt = O.OOlto, to being the natural unit 
time introduced in items (1) and (2) above. After a short 
transient (^ \{)^dt), still much bigger than the decay time 
To of self-correlations (that appears to be tq = lO'^dt), 
the system reached stationarity, in the sense that the 
instantaneous values of observables (e.g. potential en- 
ergy, Lyapunov exponents) agree with the correspond- 
ing asymptotic values within the statistical error of the 
asymptotic values themselves. After this transient we 
started recording values pi, i = 1,. . . ,J\f, of the variable 
p{x) (defined in Eq. 6), integrating the entropy produc- 
tion rate (Eq. 37) on adjacent segments of trajectory of 
length To = lOOdt = O.lto. Note that the length of the 
time interval over which we averaged the entropy produc- 
tion rate was chosen as equal to the mixing time, consis- 
tently with the discussion in Remark (4) of sec. II D. 

In conclusion, from each simulation run at fixed T and 
E we obtain Af ~ 10"" values pi of p{x) which are the 
starting point of our data analysis. The value of cr+ is 
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estimated by averaging the entropy production rate over 
the whole trajectory. 

From a shorter simulation run we measured also the 
Lyapunov exponents of the map S using the standard 
algorithm of Bcncttin et al. [25, 28]. 



D. Remarks 

To conclude this section, we note that the WCA poten- 
tial has a discontinuity in the second derivative. Thus, 
one should be concerned with the possibility that the er- 
ror of our discretization is not 0{dt^) over the qi's on a 
single time step, as it should be for potentials V ^ C*. To 
check that this is not the case (or that at least this does 
not affect our results) we made two independent tests: 

1. we simulated a system similar to model I but with 
a potential V € C'* and we obtained qualitatively 
the same results; 

2. we simulated model I using an adaptive step size 
algorithm [24]; this kind of algorithms adapt the 
step size dt during the simulation in order to keep 
constant the difference between a single step of size 
dt and two steps of size dt/2. If the precision of 
our discretization changed at the singular points of 
the potential, the time step should change abruptly 
during the simulation, while we observed a practi- 
cally constant time step during the simulation. 

Hence, we have evidence of the fact that the (isolated) 
singularities of the potentials do not produce relevant 
effects on our observations; this is probably due to the 
fact that the set of singular points of the total potential 
energy V(q) has zero measure w.r.t. the SRB measure. 



IV. DATA ANALYSIS 

In this section we will discuss in detail the procedure we 
followed to analyze the numerical data. As an example, 
we will discuss the data obtained from the simulation of 
model I a,t E = 5. As discussed in the previous section, 
from the simulation run we obtain a set Vo = {pi}i=i...j\f 
of values of the variable p{x) that correspond to r = tq 
and are measured on adjacent segments of trajectory. We 
recall again that tq = 0.1 = lOOdt is of the order of the 
mixing time, i. e. the time scale over which the correlation 
fimctions {e.g. of density fluctuations) decay to zero. 

Probability distribution function - From the dataset Vo 
we construct the histograms tTt(j>) for different values 
of T = utq as follows: the values of pix) for r = utq are 
obtained by averaging n subsequent entries of the dataset 
Vo, we obtain a new dataset V„ = {p\"''}j=i...j\f/n such 
that »^"^ = »i. Finally, from the dataset 

Vn the histogram of t^tIp) is constructed for r — nro; the 
errors are estimated as the square roots of the number of 
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FIG. 1: Model I at S = 5: the function gr{p) = Coc(p) -I- 
0((p— l)^ /t) for different values of r. 




FIG. 2: Model I at i5 = 5: the maximum Pt of Ct{p) as 
a function of 1/r. The full line is the prediction of Eq. 30, 
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counts in each bin. The function Cr(p) is then defined as 

Crip) = T-^\0gTTr{p). 

Shifting of the maximum - By fitting the function Ct- (p) in 
p € [—1,3] with a sixth-order polynomial we determine 
the position of the maximum pr within an error that, 
since Sp is the length of a bin, we estimate to be 5p/2. 
Then, we construct the function rjrip) = Ct(p — 1 + Pt) 
(see Eq. 32) which is expected to approximate the limit- 
ing function Coo(p) with error 0{(p — l)^/r). The func- 
tions rir{p) are reported in Fig. 1 for different values of r. 
We observe a very good convergence for r > 5.0 = SOtq. 

By a fourth-order fit of the so-obtained limiting func- 
tion Coo(?') around p = 1 we extract the coefficients 



9 




FIG. 3: Model I at _E = 5: the estimate of the function Coo(p) 
(open circles). In the same plot Coo(— p) +po"+ (filled squares) 
is reported. In the inset, the interval p G [—2, 2] where the 
data overlap is magnified. The full line is the Gaussian ap- 
proximation, \c}x} {p — 1)^. The plot shows that the Gaussian 
is not a good approximation in the interval [—2, 2]. The valid- 
ity of the Fluctuation Relation in the same interval is shown 
by the overlap of the open circles and filled squares. 



(j^ = -0.287 and C^' = 0.149 in order to test the cor- 
rectness of Eq. 30. In Fig. 2 we report p^- The full line is 
the prediction of Eq. 30, that is indeed verified for r 10. 
This result confirms the analysis of section II. 

Graphical verification of the fluctuation relation - From 
the previous analysis we can conclude that the function 
ryr(p) for T = 5.0 provides a good estimate of the func- 
tion Coo(p) for p e [^2,4] (see Fig. 1); thus, we can use 
this function to test the fluctuation relation, Eq. 4, in 
this range of p. In Fig. 3 we report the estimated func- 
tions Coo(p) and Coo(~p) An excellent agreement 
between the two functions is observed in the interval 
p e [—2,2] where our data allows the computation of 
both Coo(p) and Coo(— p)- Note that in this range oip the 
function Coo(p) is not Gaussian, see the inset of Fig. 3. 

Quantitative verification of the fluctuation relation - The 
translation of the function Qt{p) is crucial to obtain a 
correct estimate of the limit C,oo (p) and to verify the fluc- 
tuation relation. In this section we will try to quantify 
this observation; as the discussion will be very technical, 
the reader who is satisfied with Fig. 3 should skip to next 
section. 

The histogram 7r„,-p (p) derived from the dataset Vn is 
constructed assigning the number of counts tTq in the a- 
th bin to the middle of the binning interval, that we call 
Pa (the latter will be an increasing function of a). The 
statistical error (JtTq, on the number of counts is y^7r^. Our 
histograms are constructed in such a way that if pa is the 
center of a bin, also —pa is the center of a bin; we call 



.(3) 



a the bin such that pa = -^Pa- There exists a value Pm 
such that for pa < Pm the number of counts in the bin a 
is smaller than m (we choose m = 4). Let us indicate by 
Pa^ the smallest value of > Pm- Hence, the histogram 
is characterized by: 

1. a bin size dp; 

2. the bin am corresponding to the minimum value of 
Pa such that the number of counts in the bin is at 
least rn; 

3. the total number M of bins such that a e [£!:„,«,„]; 
for these values of pa, both 7r^(p) and T:r{—p) can 
be computed and they can be used to verify the 
fluctuation relation. 

The function Cr(p), derived from the histogram, is speci- 
fied by a set of values (pa, Ca,S(^a) for each bin a, where 
= logTTQ and the error SC,a has been defined by 



l&Ka 
T TTa 



(44) 



A quantitative verification of Eq. 4 is possible defining 
the following function: 



X = T7 



1 

M 



E 

a — Otrr 



(45) 



The value of x is the average difference between (^^ (p) and 
Ct{—p) +P<^+ in units of the statistical error. Translating 
p of a quantity a5p/2, a G Z, corresponds to shifting 
the histogram, i.e. to consider a new histogram [pa + 
a5p/2,C,a,5C,a)- This preserves the property that \i Pa is 
the center of a bin, also —pa is the center of a bin; we 
call a{a) the new value of a such that Pa{a) + a5p/2 = 
~{Pa + aSp/2). Also, the number Ma of bins such that 
a{a) £ [am,am{a)] depends on a. We define 



1 



am(a) 

E 

a—a.m 



{Ca - Ca(a) ~ (Pa + aSp/2)a+) 



N 2 



(46) 

We shall use the criterion that the fluctuation relation is 
satisfied if x 5; 3, which means that Coo(p) and Cooi~p) + 
p(T+ differ, on average, by less than 3 times the statistical 

error 

V {Kip)) + (SCi-p)) ■ The function x(a) for the 
case of model I at i? = 5 is reported in Fig. 4. 

The minimum of x is assumed between a* — 1 and 
a* -|- 1 = 2 and an upper limit for the value of x at the 
minimum is x(l) = 3.5. We estimate the translation that 
minimizes x as (5o = (a* + 0.5)Sp/2 = 1.5 • 0.093 = 0.140, 
and to this estimate we attribute an error ±Sp/2, where 
Sp — 0.186 is the size of a bin. On the other hand, we 
have seen above that, in order to shift the maximum of 
i^r (p) in p = 1 , one has to translate p by a quantity 6 = 
1 — p = 0.215. The consistency of our analysis requires 
that 6 and So coincide within their errors, i.e. that the 
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FIG. 4: Model I at S = 5: the function x(a)- The fuU 
hne corresponds to x = 3. The arrow indicates the interval 
Jo ± Jp/2 (note that its length is 2 in units of a) into which 
the minimum of x can be located within the accuracy of the 
histogram. 
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FIG. 5: Model I at -B = 5: the function X(a). The horizontal 
arrow marks the interval where the minimum of x is located, 
see Fig. 4. The vertical arrow indicates the error &X on the 
value X = 1 which is estimated as &X = 2(X(2) The 
optimal slope of the fluctuation relation without the transla- 
tion would have been X(a = 0) ~ 0.85. 



intervals J±5p/2 and (5o±(5p/2 overlap, or in other words 
|(5 — (5o| < In the present case 0.075 = |(5 — (5o| < Sp — 
0.186, then S and Sq coincide within the errors. This 
means that the translation of p brings the maximum of 
Ct{p) in p = 1 and, at the same time, minimizes the 
dilference between rir{p) and rjri—p) + p<J+, where rjr 
is our finite time estimate of Coo(p)- The value x(^*) 
quantifies this difference and is a first estimate of the 
precision of our analysis. 

Another estimate of the precision of our analysis can 
be obtained as follows. We define a parameter X as the 
slope of Coo(p) — Coo(~p) as a function of p<tj^: 

Coo{p) ^ Coo{-P) + XpG+ (47) 

The fluctuation theorem predicts X — 1, but other values 
of X are possible under different hypothesis, see [7, 10, 



29, 30]. We define a function x^(a, -'^) as 

(48) 

and for each value of a we calculate the optimal value 
of X, X(a\ by minimizing ^{a, X\ The function X{a) 
is reported in Fig. 5. As the shift of the maximum b is 
between a = 1 and a = 2, we see that the slope X is 
compatible with one. Moreover, as the natural error on 
p is the size of a bin 8p, we assign to the value A" = 1 
a statistical error 8X = 2(A(2) - A(l)) = 0.22. Note 
again that without the translation of p the optimal slope 
would be A ~ 0.85, incompatible with Eq. 4. 

Discussion - From the present analysis, we can conclude 
that: 

1. the translation shifting the maximum of Ct(p) to 
p = 1 at the same time minimizes the difference 
between ririp) and rir{—p) +pa+, where rjr is our 
finite time estimate of Coo! this proves the consis- 
tency of our theory of finite time corrections; 

2. without the translation of p (that corresponds to 
a = 0), the function Ct(p) for r ^ 5.0 rfo not sat- 
isfy the fluctuation relation, as x(a = 0) = 11 and 
A(a = 0) = 0.85; 

3. the function r]r{p) = Ct{p ~ b) satisfies the fluctu- 
ation relation with x ~ 3 and an error of about 
20% on the slope A: both quantities measure the 
accuracy of our analysis. 

Thus, the check of the fluctuation relation relies crucially 
on the translation of the function (^r (p) that has been dis- 
cussed in section II. By considering larger values of r one 
could avoid this problem (as S ~ t"^); however, as one 
can see from Fig. 1, for r > 5.0 the negative tails of Cr (p) 
are not accessible to our computational resources. The 
computation of the finite time corrections is mandatory 
if one aims to test the fluctuation relation at high values 
of the external driving force. 

Summary of the data analysis - To conclude, we sum- 
marize the procedure we follow to analyze the data of a 
given simulation run: 

1. we determine a value of r such that Cr(p) appear 
to be close to the asymptotic limit Coo(p); 

2. we determine the maximum p of Ct{p) by a sixth- 
order polynomial fit around p = 1, in an interval 
as big as possible compatibly with the request that 
the x^ from the fit is less than ^ 10; 

3. we shift the histogram of an integer multiple a of 
the half bin size Sp/2 and compute the function 
x(a) according to Eq. 46. We determine the value 
a* such that the minimum of x(fl) is assumed in the 
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FIG. 6: Model I: mobility fi as a function of the driving force 
E. The full line is the equilibrium diffusion coefficient D 
divided by the temperature. Deviations from the linear re- 
sponse are observed around E = 5. The error bars are of the 
order of the dimension of the symbols. Studying fi{E) for val- 
ues of E bigger than those shown in the figure, one can verify 
that the mobility increases up to a value jj,max , reached in cor- 
respondence of _E ~ 45. For values of E bigger than E ^ 45, 
the mobility begins to decrease essent ially followin g the lim- 
iting curve TJt/{NE), where Jt = \/T{d~ 1/N)N/T is the 
maximum allowed value of the current (saturation value). 
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FIG. 7: Model I: Lyapunov exponents for E = 5 (top) and for 
E = 25 (bottom). For each panel, the upper and lower dots 



are the two paired exponents A 

dot is their average (A^^' + A^ ' 
the dashed line is at A = 0. 



and Xj 
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and the middle 
The full line is a+/2Nd, 



interval [a* , a* + I]: the consistency of our analysis 
requires that 5 = 1 — p and 6o = (a* -t- Q.5)5p/2 
coincide within their errors (i.e. |(5 — < Sp); 

4. The value x* = ™-^^[xi'^*)iX{(^* + 1)] is an upper 
limit for the value of x at the minimum. The num- 
ber of bins min{Ma*, Ma'+i} involved in this esti- 
mate will be called M*; 

5. we compute the error 6X — 2{X {a* + 1) — X (a*)) . 

The relevant quantities r, S, Sq, \5 — do\, 5p, M*, x* and 
SX for model I are reported in table I for different values 
of the external force E. 



V. NUMERICAL SIMULATION OF MODEL I 

We will now discuss systematically the numerical data 
obtained from the simulation of model I (defined in sec- 
tion III) at different values of the driving force E. In 
Fig. 6 we report the mobility ^i{E) — T{J)e/{NE), i.e. 
the l.h.s. of Eq. 11 times T/N , as a function of E. The 
current J(p, q) has been defined in Eq. 38. From the 
Green-Kubo relation, Eq. 8, we have [9] 

(49) 



where D is the equilibrium diffusion coefficient, 

D = /nn ^ Y.^m) - q^{mE=^ (50) 

i 

Deviations from the linear response are observed and 
H{E) - D/T + 0{E^) above E ^ 5. 

In table I we report the main parameters that result 
from the data analysis (as discussed in the previous sec- 
tion) for some selected values of E. The value |(5 — (5o| 
is always less than Sp, consistently with our discussion 
above, except for E = 12.5 where, however, the relative 
difference between the two quantities is small {^-^ 9%). It 
can be noted that 6 is systematically bigger than 6q . This 
could be due to the fact that the error terms 0((p— l)^/r) 
or o(l/r) that we are discarding likely produce a system- 
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T 
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\5~So\ 


5p 


M* 


X* 


SX 


2.5 


5.0 


0.194 


0.272 0.183 


0.089 


0.244 


43 


2.2 


0.24 


5.0 


5.0 0.810 


0.215 0.139 


0.076 


0.187 


20 


3.5 


0.22 


7.5 


4.0 


1.945 


0.197 0.116 


0.081 


0.116 


18 


2.8 


0.18 


10.0 


2.5 


4.044 


0.262 0.151 


0.111 


0.122 


17 


4.4 


0.20 


12.5 


2.5 


7.090 


0.257 0.137 


0.120 


0.111 


8 


3.5 


0.28 



TABLE L Model I: results of the data analysis for some se- 
lected values of E. All the quantities are defined in section IV. 
For E > 12.5 the negative tails of the distribution are not ac- 
cessible to our numerical simulation. 
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atic shift in d or in Sq ; or that the velocity of convergence 
of (p) is not the same on the negative or on the positive 
side (because numerically is much more difficult to ob- 
serve big negative fluctuations of a than the positive ones 
- and the Fluctuation Relation provides a quantitative 
estimate of the relative probabilities). At the moment, 
because of the level of precision of our simulations, we 
are not able to investigate this problem in more detail, 
see also Remark (3) in section II D. On increasing the 
value of £', we are forced to decrease the value of r we 
use for the analysis as, for longer r, the negative tail of 
the distribution Ct{p) becomes unobservable. This can 
be seen as the number M* of bins used for the compu- 
tation of X decrease on increasing E; above E = 12.5 
it is impossible to find a value of r such that (t{p) is 
close to the asymptotic limit and the negative tail is ob- 
servable. Thus, the fluctuation relation cannot be tested 
above E — 12.5 with our computational power. However, 
we are able to check the fiuctuation relation in the region 
E > 5 where deviations from the linear response are ob- 
served. Moreover, the estimated distributions Coo(p) are 
very similar to the one reported in Fig. 3: in particular, 
they are not Gaussian in the investigated interval of p 
(also for i? < 5, in the linear response regime). 

Finally, in Fig. 7 we report the measured Lyapunov 
exponents of the model ior E = 5 and = 25. For this 
system, the Lyapunov exponents are known to be paired 
[25, 31, 32] like in Hamiltonian systems and the average 
of each pair is a constant equal to a+/2Nd. For E — 5, 
each pair is composed of a negative and a positive expo- 
nent. This means that the attractive set is dense in phase 
space [10, 30] and the chaotic hypothesis is expected to 
apply to the system yielding a slope X = 1 in the fluctu- 
ation relation, as confirmed by our numerical data. The 
same happens up to £^ ^ 20. Above E = 20, there is a 
number D of pairs composed by two negative exponents 
(for _E = 25 we get D = A, see Fig. 7). In this situation, 
the slope X in the fluctuation relation is expected to be 
given by X = 1 - D/Nd [29, 30]. Thus, for E = 25 one 
expects X ^ 0.75. Unfortunately, as discussed above, 
above E = 12.5 we did not observe negative fluctuations 
of the entropy production, and this prediction could not 
be tested in our simulation. 



VI. NUMERICAL SIMULATION OF MODEL II 

Model II differs from model I in the dimension d = 3, in 
the larger number of particles iV = 20, and because it is a 
binary mixture of two types of particles. Binary mixtures 
are frequently used as models for numerical simulations 
of supercooled liquids as they avoid crystallization also 
at very low temperature on the "physical" time scales 
(i.e. on the time scales of numerical experiments); for 
these systems, at low temperature deviations from the 
linear response are observed also for very low values of 
the external driving force. 

In Fig. 8 we report the equilibrium diffusion coeffi- 
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FIG. 8: Mobility as a function of the temperature T and of 
the driving force E for model II. The circles correspond to the 
equilibrium diffusion coefficient divided by the temperature. 
Deviations from the linear response are observed for E > 3; 
they become larger on lowering the temperature, as D — > 0. 




FIG. 9: The estimate of the function Coo(p) (open circles) for 
model II with T = 1.1 and E — 3. In the same plot Coo(— p) + 
pa+ (filled squares) is reported. In the inset, the interval 
p G [—1.5, 1.5] where the data overlap is magnified. The full 

line is the Gaussian approximation, Coo(p) = |Coo''(p — 1)^- 
The data have been obtained from the histogram of nrip) 
with T = 2.5 (see table II). 



cient D (divided by the temperature T) and the mobility 
(for different values of E) as functions of the tempera- 
ture. Even though the number of particles is very small, 
on lowering the temperature the systems approaches the 
supercooled state and D becomes very small around 
T ^ 0.5. Slightly above this temperature, i.e. around 
T = 1, strong deviations from the linear response are 
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observed for E > 3, where the entropy production (7+ is 
still close to 0. Some values of (7+ are reported in table II; 
to compare these values with those obtained for model I 
one should note that a+ is an extensive quantity. Thus, 
the entropy production per degree of freedom, (T_|_ /2Nd, 
is much smaller in model II than in model I. 

In table II the results of the data analysis outlined in 
section IV are reported. For i? < 6 we obtain a very 
good agreement of the data with the predictions of the 
fluctuation relation and with the theory of finite time 
corrections discussed in section II. For = 10 it is very 
difficult to observe negative fluctuations of p with our 
computational power; see e.g. the result of the analysis 
for £■ = 10 and T = 1.9, where only M* = 7 bins where 
available and we were forced to use r = 0.2, of the order 
of the mixing time Tq. In Fig. 9 we report the estimated 
function Coo(p) obtained for T = 1.1 and E = 3 from the 
data with r = 2.5. Strong deviations from the Gaussian 
behavior are observed in the accessible range of p (see 
the inset of Fig. 9). A similar behavior of Coo(p) is ob- 
served in correspondence of all the values of E and T we 
investigated (those listed in Table II) : in particular in all 
these cases highly non Gaussian behaviors are observed 
in the accessible range of p. 

The Lyapunov spectrum for this system is very similar 
to the one reported in the upper panel of Fig 7. Pairs of 
two negative exponents were observed only for E = 10 
at T < 1.3, where, as in the case of model I, a+ is too 
large to allow for a verification of the modified fluctuation 
relation expected in this case, see the discussion at the 
end of section V. 



VII. CONCLUSIONS 

We tested the fluctuation relation, in our opinion quite 
successfully, in a numerical simulation of two models of 



T E 


T a+ 5 So \S~do\ 6p M* x* 5X 


0.9 1 
0.9 3 


3.0 0.209 0.453 0.334 0.119 0.223 68 1.9 0.19 
3.0 2.615 0.286 0.264 0.024 0.132 15 1.0 0.23 


1.1 1 
1.1 3 
1.1 6 


4.0 0.233 0.231 0.126 0.105 0.126 79 1.7 0.24 
2.5 2.493 0.217 0.238 0.021 0.087 30 1.0 0.12 

1.5 13.32 0.113 0.230 0.117 0.092 7 1.1 0.21 


1.5 1 
1.5 3 
1.5 6 


3.0 0.230 0.179 0.140 0.039 0.140 86 0.9 0.13 
2.5 2.227 0.145 0.123 0.022 0.082 33 4.7 0.18 
0.5 52.14 0.074 0.130 0.056 0.052 11 0.6 0.10 


1.7 1 


3.0 0.221 0.127 0.141 0.014 0.283 49 1.0 0.26 


1.9 3 
1.9 6 

1.9 10 


2.5 1.981 0.106 0.122 0.016 0.122 26 0.8 0.12 
0.4 43.52 0.078 0.126 0.048 0.085 14 1.7 0.11 

0.2 139.0 0.079 0.135 0.056 0.039 7 0.8 0.10 


2.1 6 


0.4 40.48 0.074 0.110 0.036 0.110 11 1.0 0.15 



TABLE II: Model II: results of the data analysis for some 
selected values of T and E. All the quantities are defined in 
section IV. 



interacting particles subjected to an external iioiiconscr- 
vative force and to a reversible mechanical thermostat. 
Our data satisfy the fluctuation relation with a x < 3 and 
an accuracy of the order of 20% also for very large val- 
ues of the driving force, where strong deviations from the 
linear response are observed, and where the large devia- 
tion function is strongly non-Gaussian. The comparison 
of our numerical data with the predictions of the fluctu- 
ation relation is done by taking into account the (lowest 
order) flnite time corrections to the distribution function 
for the fluctuations of the phase space contraction rate. 
This is crucial: if we did not take into account such cor- 
rections the fluctuation relation would be violated within 
the precision of our experiment. 

In order to compute the finite time corrections, we 
proposed an algorithm which allows to reconstruct the 
asymptotic distribution function from measurable quan- 
tities at finite time, within a given precision. Our theory 
of the corrections relies on the symbolic representation 
of the chaotic dynamics, therefore it is applicable if one 
accepts the Chaotic Hypothesis. 

Our interpretation of the numerical results is that the 
chaotic hypothesis can be applied to these systems, also 
very far from equilibrium, and in particular the fluctu- 
ation relation is satisfied even in regions where its pre- 
dictions measurably differ from those of linear response 
theory. 

Our theory of finite time corrections for the analysis 
of our numerical data could in principle be of interest for 
real experimental settings where non Gaussian fiuctua- 
tions for the entropy production rate are observed, see 
[35, 36]. 

However it should be stressed that in a real experiment 

there arc some technical differences with respect to our 
numerical simulation which could in some cases make in- 
applicable our analysis, namely: 

(i) usually the noise in the large deviation function for 
the entropy production rate in a real experiment is much 
bigger than in a numerical experiment, and it is likely 
that the translation in Eq. 32 computed as the ratio 
C^^V(C^^^)^ is not measurable within an error of some 
percent; 

(a) usually in a real experiment the accessible time scales 
are naturally much bigger than the microscopic ones so 
that, if the negative fluctuations of the entropy produc- 
tion rate arc observable at all, one is automatically in 
the asymptotic regime, where the finite time corrections 
should be negligible; 

(Hi) a usual problem in a realistic setting is that there 
is no clear connection between the "natural" thermo- 
dynamic entropy production rate s = W/T [W is the 
work of the dissipative external forces and T is the tem- 
perature) and the microscopic phase space contraction 
rate, for which a slope AT = 1 in the fluctuation relation 
C,{p) — C,{—p) = Xa+p is expected; so, often one measures 
an X 7^ 1 and correspondingly one defines an effective 
temperature Te// — T/X giving a natural connection 
between the effective thermodynamic entropy production 
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rate Seff — W/T^/ / and the phase space contraction rate, 
see [7, 35, 36]; in such a situation (where an adjustable 
parameter X appears) it makes no sense to apply our 
analysis, which is sensible only if one wants to compare 
the experimental data with a sharp prediction about the 
slope X in the fluctuation relation. 

A big open problem we are left with is trying to under- 
stand how the fluctuation relation is modified for values 
of the driving force so high that the attractive set is no 
longer dense in phase space. It is expected, [29], that in 
such a case Coo (p) ~ Coo i^p) is still linear, but the slope is 
X(T_|_, with X given by the ratio of the dimension of the 
attractive set and of that of the whole phase space. An 
estimate of such quantity can be given via the number 
of negative pairs of exponents in the Lyapunov spectrum 
[29, 30]. Unfortunately negative pairs begin to appear in 
the Lyapunov spectrum only for values of the external 
force so high that no negative fluctuations are observable 
anymore. We hope that future work will address this 
point. 



APPENDIX A: A LIMIT THEOREM 



let af, TT^ and 11-^(6; I) be defined as above. Then, for 
a sufficiently sm,all e > 0, there exists an analytic "rate 
function" Ct(p) such that 



Ur{e;I) 



lim 



1 . 



(Al) 



Ct{p) is defined by: 



Ct{p) + - log 



sinh[eAp/ {2a -^ 



eK/{2a+) 



C(p) 



Crip) 



-Ze(A;) + A>a+ - — log[— ( - 



(A2) 



and Ap is the inverse of p{X) = z^(A)/(t+. The function 
C^ ip) has the following property: if A C I is an interval 
of size around a point p/^, then: 



lim ^-(^i^) 



(A3) 



In this section we prove Eq. 23-24. We reproduce in 
detail the proof in the case p is the average of indepen- 
dently distributed discrete variables erf, assuming values 
in eZ, for some small mesh parameter e; then we discuss 
how this can be applied and adapted to the situation 
considered in section IIC and subsequent sections. 

Let us introduce some definitions. Let (Ti, i € N, be 
independent continuous random variables with identical 
distributions TT{dai) with positive variance Sa^ > 0, sup- 
ported on the finite interval [s_,s+]. Let us assume 
that TT{d(Ti) gives positive probability to any finite in- 
terval contained in [,s_,s_|_]. Let ■n\{dij) be the weighted 
distribution ■n\{d<7) — e~^'^7r((iCT)/ / e~'^'^7r(dcr) and let 
us define 2:oo(A) = — log/ e~'^'^7r(dcr) and cr+ = 2:^(0). 
Note that the assumption that ■n{dai) gives positive prob- 
ability to an interval of a in [s_, s+\ implies that for any 
finite A also 'K\{d(j) has positive variance —z'^{X) > 0. 

Also, given e > (with the property that s+ — .s_ = 
NgS for some integer N^), let us consider the discretiza- 
tion of ai on scale e, call it erf: erf will be a discrete 

variable assuming the values sl'^=s^ + (fc — i)e, k = 
1, . . . , Ng, with probabilities 7r^{s%) = Prob(crf = s%) = 
f e , e n(da). The assumption that iridai) gives positive 

probability to any finite interval contained in [s_, s+] im- 
plies that T^^{s\) > for any e and k. Let also Zs{X) = 
-logE^ie-'^^^Tr^sl) iindnl{s%)=n%s%)e-^^l+^^W. 
Note that, since 7r^(s|) > for any k, for any finite A one 
has -4'(A) > 0. 

If Pr = 7^ Xli=i ^'^d ni-(e;/) is the probabihty 
that p^ belongs to the finite interval I, the following 
theorem holds. 

Theorem: Given a finite interval I C (s_,s+), 



Proof Let us introduce the auxiliary variable q = 
7^Sl=i'7i> where rji are i.i.d. discrete random vari- 
ables, with distribution 7r|(s|). Let us call II^{e;qo) 
the probability that q assumes the value qo (z I, with 
qo(y+ = sl/r for some A; e N, and note that n°(e;go) is 
identical to the probability that Pr = qo- By definition 
Il^{qo) and Il^{qo) are related by: 



n^(e;ao) 



e-^"""'"n';(r:ryo) 
[E.e-^^W^(4)]^ 



(A4) 



Now, a local form of central limit theorem (Gnedenko's 
theorem, see pag. 211 of [33]) tells us that, if g is localized 
near its mean value, that is if \qa-^- — z'g{X)\ < ^ for some 
finite M, then JI^{e;qo) is asymptotically equivalent to 
the Gaussian with mean z'^{X) and variance —z'J{X), in 
the sense that 



n^(9o) 



V27rr(-4'(A)) 



^(1 + 0(1)) 



(A5) 

for any go s.t. |ga+-z^(A)| < ^ [37]. 

So, given A^^ s.t. ^(^90) = 1"'^+ (^^^^ Ko exists, is 
unique and is an analytic function of qQ, by the remark 
that —2" (A) > for any finite A and Ze{X) is an analytic 
function of A), using Eq. A5 we see that Eq. A4 can be 
restated as: 



n?(e;go) 



^2nT{-z'J{XlJ) 



eKoi°''+-'-^^(Ko\l + o{l)) 
(A6) 



Now, by the definition of Ct(p) i^i ^1- ^2, we see that 
the r.h.s. of the last equation is equal to :p^e^^?^'°^(l -|- 
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o(l)). Finally, the statement of the Theorem follows by 
the remark that 

rCApo) = dpe"^-(P)(l+o(l)) . (A7) 



rf7+ 



In fact the integral in the r.h.s. of the last equation is 
given by 



_crC.fon) 2sinh[C;(po)£/(2a+)] / C{po)e^ \ 

<(P0) ^ ^ T V 

(A8) 

and in the last expression one has to note that 

Cipo) = [CrYiPo) + o(i) = a;„ + o(i). ■ 

A first Remark to be done about the Theorem above is 
that, in order to define a "universal" rate function in 
terms of quantities depending only on 2;oo(-^) (instead 
of quantities depending on the "non universal" function 
^^(A), which explicitly depends on the discretization 
step e), it would desirable to perform (in a sense to be 
precised) the continuum limit e ^ 0. To this regard, 
we can note that the only point where in the proof 
above we really used the fact that e is a constant (i.e. 
is independent of t) was in using Gnedenko's Theorem, 
sec [33]. However, by a critical analysis of the proof 
of Gnedenko's Theorem, one can realize that it is even 
possible to let e = e,- go to with r; the velocity with 
which Er is allowed to go to depends on the details 
of the distribution ^{da). So we can even study the 
probability distribution of Pr on a scale ~ er/r: if we 
introduce bins A,- of size 0{er/T) and we define liri^r) 
to be the probability that Pr = J2i belongs to the 
bin Aj- centered in po, we can repeat the proof above to 
conclude that 



lim 



n.(A.) 



1 



(A9) 



where Cr satisfies the equation: 



(Ap) = -zU\) + Appa+ - ^ log[- ( - f^)] 

It TV cr_|_ / 

(AlO) 



and Ap is the inverse of p(A) = z^(A)/ct+. 

Another point to be discussed is that in the Theorem 
above we assumed the cji to be independent. This is not 
the case for the variables cr(S'*-) of sec. II. However, if, 
as discussed in Remark (3) of sec. II D, we choose the 
time unit to be of the order of the mixing time, the 
variables cr(S'*-) have (by construction) a decorrelation 
time equal to 1, and the analysis of previous theorem 
can be repeated step by step in order to construct the 
probability distribution oi p = -^-^ (7(5'* •). The only 
differences are that: (1) t2:oo(A) should be replaced 
by TZr{X) = -log/e-^J""+'^n^(dp) throughout the 
discussion; (2) instead of Gnedenko's theorem one has 
to apply a generalization of Gnedenko's to short ranged 
Gibbs processes, to be proven via standard cluster 
expansion techniques (see for instance [34] for a proof of 
a generalization of Gnedenko's theorem to a short ranged 
Gibbs process in the context of non critical fluctua- 
tions of the phase separation line in the 2D Ising model). 

The conclusion is that, if the bins A in sec. II C are 
chosen of size £t/t, the probability of the bin A centered 
in Pa is asymptotically given by Tr{p e A) ~ e'^^^^P^^ (in 
the sense of Eq. 23) and Ct{pa) can be interpolated by 
an analytic function of p that in fact satisfies Eq. 24. 
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